

***Table1

use "$working_data/Air_Data.dta", clear
sum2docx manipulate_pm25 did wind humidity temper lnrain using "$tables/result.docx", ///
    append stats(N mean(%6.2f) sd(%6.2f) min(%6.2f) max(%6.2f)) ///
    title("")
	
***Table 2
use "$working_data/Air_Data.dta", clear
global controls "wind humidity temper lnrain lngdp lnpergdp second fiscal"

global fe "absorb(citycode yrmt)"
global vc "vce(cluster citycode)"

reghdfe manipulate_pm25 did,noabsorb $vc
reghdfe manipulate_pm25 did $controls, noabsorb $vc
reghdfe manipulate_pm25 did $controls, absorb(citycode) $vc
reghdfe manipulate_pm25 did $controls, $fe $vc

***Table3
use "$working_data/Air_Data.dta", clear
global controls "wind humidity temper lnrain lngdp lnpergdp second fiscal"
global fe "absorb(citycode yrmt)"
global vc "vce(cluster citycode)"

reghdfe manipulate_pm25 c.did##c.target $controls, $fe $vc
reghdfe manipulate_pm25 c.did##c.target $controls, $fe $vc
reghdfe manipulate_pm25 c.did##c.haze_baidu $controls, $fe $vc
reghdfe manipulate_pm25 c.did##c.haze_baidu $controls, $fe $vc 
 
***Table4
use "$working_data/Air_Data.dta", clear
merge m:1 city  using "$temp_data/policy_city.dta", keep(3) nogen
gen reform     = (yrmt >= policy)
global controls "wind humidity temper lnrain lngdp lnpergdp second fiscal"
global fe "absorb(citycode year)"
global vc "vce(cluster citycode)"

reghdfe manipulate_pm25 reform, noabsorb $vc
reghdfe manipulate_pm25 reform $controls, noabsorb $vc
reghdfe manipulate_pm25 reform $controls, absorb(citycode)
reghdfe manipulate_pm25 reform $controls,$fe $vc 
  
***Table5
use "$working_data/Air_Data.dta", clear
gen month = month(dofm(yrmt))
inlist2 city, values(北京市, 天津市, 上海市, 重庆市) name(city4)
reghdfe manipulate_so2 did $controls, $fe $vc
reghdfe manipulate_pm25 did $controls if mi(city4), $fe $vc
reghdfe manipulate_pm25 did $controls if year != 2020, $fe $vc
bootstrap, reps(100): reghdfe manipulate_pm25 did $controls, $fe
reghdfe manipulate_pm25 did  inspect_did  $controls , $fe $vc

		
***Table6
use "$working_data/Air_Data.dta", clear
global controls "wind humidity temper lnrain lngdp lnpergdp second fiscal"
global fe "absorb(citycode yrmt)"
global vc "vce(cluster citycode)"

reghdfe lncases  did $controls ,$fe $vc
reghdfe  manipulate_pm25  lncases  $controls  ,$fe $vc

reghdfe lnemployee did $controls, $fe $vc 
reghdfe manipulate_pm25 lnemployee $controls, $fe $vc



***Figure 1
use "$working_data/intro.dta", clear
coefplot
    (m1, label("PM2.5(satelite)") msize(*.5) m(o) mcolor(gs0) ciopt(lc(gs0)))
    (m2, label("PM2.5(official)") msize(*.5) m(T) mcolor(gs9) ciopt(lc(gs9))), 
    drop(_cons)
    aseq vertical
    ytitle("Coefficient", size(*.8))
    xlabel(0(12)72 0 "2014-01" 12 "2015-01" 24 "2016-01" 36 "2017-01" 
           48 "2018-01" 60 "2019-01" 72 "2020-01", labsize(*.85) nogrid)
    ylabel(-60(10)20, angle(0) labsize(*.85) nogrid)
    xline(0(3)81, lstyle(grid) lwidth(thin) lcolor("235 235 235"))
    yline(-60(5)20, lstyle(grid) lwidth(thin) lcolor("235 235 235"))
    yline(0, lp(dash))
    legend(row(1) ring(1) size(*.9) position(6))
    graphregion(fcolor(gs16) lcolor(gs16)) 
    plotregion(lcolor(gs16) lwidth(*0.9))
    scheme(s2mono);
#de cr
graph export "$figures/preface.svg", width(1000) replace



***Figure 3
use "$working_data/Air_Data.dta", clear
duplicates drop provcode yrmt, force
xtset provcode yrmt

panelview manipulate_pm25 did, i(provcode) t(yrmt) type(treat) bytiming
    xtitle("Date", size(*.5))
    ytitle("Province", size(*.5))
    xlabel(1 "2014-2" 11 "2014-12" 23 "2015-12" 35 "2016-12" 47 "2017-12" 59 "2018-12" 71 "2019-12" 83 "2020-12")
    ylabel(1 "Chongqing" 2 "Hebei" 3 "Jiangsu" 4 "Shandong" 5 "Qinghai" 6 "Hubei" 7 "Fujian" 
           8 "Jiangxi" 9 "Tianjin" 10 "Shanghai" 11 "Shaanxi" 12 "Guangxi" 13 "Guangdong" 14 "Inner Mongolia"
           15 "Ningxia" 16 "Jilin" 17 "Henan" 18 "Guizhou" 19 "Yunnan" 20 "Anhui" 21 "Sichuan"
           22 "Gansu" 23 "Zhejiang" 24 "Beijing" 25 "Liaoning" 26 "Hunan" 27 "Shanxi" 28 "Heilongjiang"
           29 "Hainan")
    mycolor(mono)
    legend(row(1) ring(1) position(6) size(*.8));
	
***Figure4a
use "$working_data/Air_Data.dta", clear
gen dist = yrmt - pilot
two (kdensity pm25_satelite if year == 2014, bwidth(12) boundary 
        color("black") lp(solid)) //  
    (kdensity pm25_satelite if year == 2015, bwidth(12) boundary 
        color("black") lp(dash))  //  
    (kdensity pm25_official if year == 2014, bwidth(12) boundary 
        color("black") lp(solid)) //  
    (kdensity pm25_official if year == 2015, bwidth(12) boundary 
        color("black") lp(dash)), //  
    xtitle("")
    xlabel(0(25)225, labsize(*.85))
    xline(0(12.5)225, lstyle(grid) lwidth(thin) lcolor("235 235 235"))
    ytitle("")
    ylabel(0(0.005)0.025, angle(0) format(%4.3f) labsize(*.85))
    yline(0(0.00125)0.025, lstyle(grid) lwidth(thin) lcolor("235 235 235"))
    legend(order(1 "satelite_2014"
              2 "satelite_2015"
              3 "official_2014"
              4 "official_2015")
           position(bottom) //  
           region(lstyle(solid) lcolor(black) fcolor(white)) // 
           row(2) col(2) size(*.8)) //  
    graphregion(fcolor(gs16) lcolor(gs16)) 
    plotregion(lcolor("white") lwidth(*0.9))
    text(0.0235 36 "Satelite " , color(black) size(*.8)) //  
    text(0.0022 150 "Official", color(black) size(*.8)); //  
    scheme(s2color) 
	
***Figure4b
two (kdensity pm25_satelite if year == 2018, bwidth(12) boundary 
        color("black") lp(solid)) //  
    (kdensity pm25_satelite if year == 2019, bwidth(12) boundary 
        color("black") lp(dash))  //  
    (kdensity pm25_official if year == 2018, bwidth(12) boundary 
        color("black") lp(solid)) //  
    (kdensity pm25_official if year == 2019, bwidth(12) boundary 
        color("black") lp(dash)), //  
    xtitle("")
    xlabel(0(25)175 190, labsize(*.85))
    xline(0(12.5)175 190, lstyle(grid) lwidth(thin) lcolor("235 235 235"))
    ytitle("")
    ylabel(0(0.005)0.025, angle(0) format(%4.3f) labsize(*.85))
    yline(0(0.00125)0.025, lstyle(grid) lwidth(thin) lcolor("235 235 235"))
    legend(order(1 "Satelite 2018"
              2 "Satelite 2019"
              3 "Official 2018"
              4 "Official 2019")
           position(bottom) 
           region(fcolor(white) lcolor(black)) //  
           row(2) 
           cols(2) //  
           size(*.8)) //  
    graphregion(fcolor(white) lcolor(black)) 
    plotregion(lcolor("white") lwidth(*0.9))
	text(0.024 36 "Satelite " , color(black) size(*.8)) //  
    text(0.0013 150 "Official", color(black) size(*.8)); //  
    scheme(s2color) 	


	
***Figure5a
#de ;
marginsplot, 
    title("")
    ciopts(color("black"))
    plotopts(
        msize(*.6)
        color("black")
        xtitle("Economic growth target", size(*.8))
        xlabel(, labsize(*.85))
        xline(5(.5)15, lstyle(grid) lwidth(thin) lcolor("235 235 235"))
        ytitle(" Marginal effect ")
        ylabel(-0.1(0.05)0.1, format(%4.2f) angle(0) labsize(*.85))
        yline(-0.1(0.0125)0.1, lstyle(grid) lwidth(thin) lcolor("235 235 235"))
        yline(0, lp(dash))
        graphregion(fcolor(gs16) lcolor(gs16)) 
        plotregion(lcolor("white") lwidth(*0.9)));
#de cr

***Figure5b
#de ;
marginsplot, 
    title("")
    ciopts(color("black"))
    plotopts(
        msize(*.6)
        color("black")
        xtitle("Baidu_index", size(*.8))
        xlabel(, format(%2.1f) labsize(*.85))
        xline(0(0.25)6, lstyle(grid) lwidth(thin) lcolor("235 235 235"))
        ytitle("Marginal effect")
        ylabel(-0.8(0.2)0, format(%4.2f) angle(0) labsize(*.85))
        yline(-0.8(0.05)0, lstyle(grid) lwidth(thin) lcolor("235 235 235"))
        yline(0, lp(dash))
        graphregion(fcolor(gs16) lcolor(gs16)) 
        plotregion(lcolor("white") lwidth(*0.9)));
#de cr 


***Figure6
 
two (scatter estimate id, c(l) color(black) msize(small))
    (scatter ul id, c(l) m(none) color(gs6) lp(dash))
    (scatter ll id, c(l) m(none) color(gs6) lp(dash)),
    xtitle("Relative time of reform implementation")
	ytitle("Regression coefficient")
    yline(0)
    xline(-1, lp(dash) lc(black))
    xlabel(-10(2)14 15, grid)
    ylabel(, grid format(%4.2f))
    legend(off)
    graphregion(color(white)) 
    plotregion(lcolor(black)); 


***Figure7
 
cap erase "$working_data/simulations.dta"
permute did beta = _b[did] se = _se[did] df = e(df_r), ///
    reps(1000) rseed(123) saving("$working_data/simulations.dta"): ///
        reghdfe manipulate_pm25 did $controls, $fe
use "$working_data/simulations.dta", clear
two (histogram beta, bin(60) finten(50))
    (kdensity beta, bw(0.0025) lp(solid)),
    xtitle("Estimator", size(*.85))
    xlabel(-0.02(0.005)0.02, format(%4.3f) labsize(*0.85))
    xline(-0.0181, lp(dash)) 
    xline(-0.02(0.0025)0.02, lstyle(grid) lwidth(thin) lcolor("235 235 235"))
    ytitle("Density", size(*.85))
    ylabel(0(30)150, angle(0) labsize(*0.85))
    yline(0(15)150, lstyle(grid) lwidth(thin) lcolor("235 235 235"))
    legend(off)
    graphregion(fcolor(gs16) lcolor(gs16)) 
    plotregion(lcolor("white") lwidth(*0.9))
    scheme(s2mono);



 











